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Abstract. We provide a Jacobian criterion that applies to arbitrary chemical reaction 
networks taken with mass-action kinetics to preclude the existence of multiple positive 
steady states within any stoichiometric class for any choice of rate constants. We are 
concerned with the characterization of injective networks, that is, networks for which the 
species formation rate function is injective in the interior of the positive orthant within 
each stoichiometric class. We show that a network is injective if and only if the de- 
terminant of the Jacobian of a certain function does not vanish. The function consists 
of components of the species formation rate function and a maximal set of indepen- 
dent conservation laws. The determinant of the function is a polynomial in the species 
concentrations and the rate constants (linear in the latter) and its coefficients are fully 
determined. The criterion also precludes the existence of degenerate steady states. Fur- 
ther, we relate injectivity of a chemical reaction network to that of the chemical reaction 
network obtained by adding outflow, or degradation, reactions for all species. 

Keywords: Jacobian criterion, multiple steady states, injectivity, stoichiometric space, 
degenerate steady state 

1. Introduction 

Multistationarity in cellular systems provides a mechanism for switching between differ- 
ent cellular responses and can be crucial for cellular decision making. Even though different 
features, such as feedback loops, are known that facilitate multistationarity in systems, it 
is in general difficult to decide whether a particular system has the capacity to exhibit mul- 
tiple steady states. Typical systems are high-dimensional and contain many parameters 
that are unknown or poorly determined. In order to determine the steady states of such 
a system, the simultaneous solutions to a large set of equations taken together with the 
unknown parameters is required. In general, this is an impractical task. Various criteria 
have therefore been developed to preclude the existence of multiple (positive) steady states. 
These criteria typically utilize the structure or qualitative features of the system [U El [l] 
or properties of the class of kinetics that are allowed [H 13 12] ■ 

It is the aim of this paper to introduce a criterion for a chemical reaction network 
taken with mass-action kinetics to preclude the existence of multiple positive steady states 
within any stoichiometry class for any choice of rate constants. The criterion is based on 
the species formation rate function and characterizes when this function is injective for 
positive concentration vectors within each stoichiometric class. If this is the case then the 
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network is said to be injective (Definition 5.1) and there cannot exist multiple positive 
steady states within a stoichiometric class. 

We provide a Jacobian criterion that characterizes injectivity for any network. The crite- 
ria is computationally tractable and extends the Jacobian criterion for fully open networks 
in [3]. In a fully open network all chemical species are assumed to flow out of the system or, 
alternatively, all species are being degraded. If the dimension of the stoichiometric space 
is less than the number of species then the Jacobian of the species formation rate function 
is always singular and the criterion for fully open networks does not apply. We replace the 
species rate formation function by a new function (Definition 4.4) obtained by replacing 
some components of the species formation rate function by equations for independent con- 
servation laws. We show that a network is injective if and only if the determinant of the 
Jacobian of the modified function does not vanish (Corollary 5.4). If this is the case then 
any positive steady state is non- degenerate (Corollary 5.5). The criterion depends only on 
the structure of the network and, therefore, is easy to implement using any computational 
algebra software. 

Any network can be seen as a subnetwork of a fully open network by adding outflow 
reactions. We show that the determinant of the Jacobian related to the original network 
can be recovered from that of the associated fully open network. Further, we show that if 
the fully open network is injective then either the original network is injective as well and 
all positive steady states are non-degenerate, or all steady states of the original network are 
degenerate (Theorem 9.1). Further, the cases for which the latter occurs are characterized 
(Corollary 8.1). 

Our work builds on previous work by Craciun and Feinberg [4J on injectivity of networks 
in the context of a continuous flow stirred tank reactor (here called fully open networks). 
In two subsequent papers Craciun and Feinberg relate these results to arbitrary networks. 
First, in |5j they show that if a fully open network does not have the capacity to admit 
multiple positive steady states, then any network resulting from removing outflow reactions 
cannot have multiple non-degenerate positive steady states within any stoichiometric class. 
Later, in [7], they provide a (sufficient) condition to ensure that in an arbitrary network 
degenerate steady states cannot exist assuming that the fully open network is injective. 
Therefore, the combination of the two results gives a criterion to preclude multiple positive 
steady states. It consists of first deciding whether the associated fully open network is 
injective and then determining if degenerate steady states can occur. However, a network 
can be injective even if the associated fully open network is not. 

The work presented here provides a direct path to preclude multiple steady states by 
avoiding the detour to fully open networks. We show that injectivity of an arbitrary network 
can be assessed regardless of the injectivity of the associated fully open network. Further, 
if the associated fully open network is injective, the occurrence of degenerate steady states 
in a network is completely characterized. 

A different route to injectivity of a fully open network was taken by Banaji et al. in 
[2] . A criterion is given that ensures that minus the Jacobian of the species formation rate 
function is a P- matrix (the definition is given in QO). It then follows from the results of 



|12j that the network is injective. Our results imply that, after changing the sign of certain 
rows, the Jacobian of the modified species formation rate function is a P-matrix and it 
follows that the network is injective as well. Using this approach and the notion of strongly 
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sign-determined matrices, Banaji et al. extend in f2] the injectivity results of Craciun and 
Feinberg for fully open networks taken with mass-action kinetics to kinetics satisfying some 
mild conditions (see also |T]). Our work is currently restricted to mass-action kinetics and 
the extension to general kinetics is currently being investigated. 

The outline of the paper is as follows. In ^we introduce some notation and the main 
definitions relating to networks and mass-action kinetics. In ^we introduce the stoichio- 
metric classes and the distinction between fully open and closed networks. We proceed in 
^to study degenerate steady states. Injectivity of networks is discussed in ^ where the 
definition of injectivity and the Jacobian criterion are introduced. In ^and ^we focus on 
open and closed networks, respectively. Section [8] provides a characterization of networks 
with only degenerate steady states. Finally, in ^we relate injectivity of open networks to 



that of closed networks and in { 10 the relationship between P-matrices and injectivity is 
discussed. We end with a few remarks including a summary (Figure [T]) of how our work 
relates to the previous work of Craciun and Feinberg. 

2. Chemical reaction networks with mass-action kinetics 

2.1. Notation. Let M+ denote the set of positive real numbers (without zero) and M+ the 
set of non-negative real numbers (with zero). Similarly, let N be the set of non- negative 
integers. Given a finite set £, let N be the semi-ring of formal sums v = YIegs ^eE, with 

A^; G N. If A^; G N for all E' € then we write u G N^. The semi-rings and are 
defined analogously. 

The ring of polynomials in £ is denoted M[f ]. The total degree of a monomial n_Bef -^"^) 
with ue a non-negative integer for all is the sum of the degrees of the variables, 
Ssef The degree of a polynomial is the maximum of the total degrees of its monomi- 
als. 

If a polynomial p vanishes for all assignments a : — )■ then p = identically. Further, 
if p is a non-zero polynomial in ]R[£^] such that the degree of each variable in each monomial 
is either 1 or zero, then all the coefficients of p are non-negative if and only \i p{a{8)) > 
for any assignment a: £ ^ M_|_ . 

2.2. Chemical reaction networks. Here we introduce the definition of a chemical re- 
action network and some related concepts. See for instance [U |9] for background and 
extended discussions. 

Definition 2.1. A chemical reaction network (or simply network) consists of three finite 
sets: 

(1) A set 5 of species. 

(2) A set C C N of complexes. 

(3) A set 7^ C C X C of reactions, such that {y,y) ^ TZ for all y £ C, and if y G C, then 
there exists y' £ C such that either {y,y') £ TZ or {y',y) € TZ. 

Following the usual convention, an element r = {y,y') £ TZ is denoted by r: y — )• y'. 
The reactant and product (complexes) of a reaction r: y ^ y' are y and y' , respectively. 
By definition, any complex is either the reactant or product of some reaction. The zero 
complex G C is allowed by definition. Reactions of the form y — ?• are called outflow 
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reactions and reactions of the form — )• y are called inflow reactions |10j . In particular, 
species inflow and species outflow reactions are reactions of the form — )• S and 5 — )• 0, 
respectively, for some S £ S. 

Let n be the cardinality of S. We fix an order in S so that S = {Si, . . . , Sn} and identify 

N with N . The species Si is identified with the i-th canonical n-tuple of N with 1 in the 
i-th position and zeroes elsewhere. Accordingly, a complex y G C is given as y = Y17=i Ui^i 
or (yi, . . . , y„). Although N is not a vector space, n-tuples f G N will be called vectors. 

Example 2.2. Consider the network with set of species S = {Si, S2, S^, Sa, S^, Sq} , set of 
complexes {Si + S-i, Si + S^, S2 + Si, S2 + S-i, S^, S^} and reactions 

Si + < ^ 5*5 *- 5i + 5*4 S2 + S4 < ^ Sq 52 + S3. 

That is, the reactions are ri : 5i + ^3 — ?■ S5, r2: S5 — )■ Si + S3, rs: S5 — )■ Si + S4, 
r4 : 5*2 + 5*4 — )• Sg , rs : Sq ^ S2 + S^ and vq: Sq ^ S2 + S^. This network is a main building 
block in protein modification systems and is known as the futile cycle. It assumes the 
Michaelis-Menten enzyme mechanism in which a substrate S3 is modified into a substrate 
54 through the formation of an intermediate complex 5*5. The reaction is catalyzed by 
an enzyme ^i. The modification can be reversed via a similar set of reactions with an 
intermediate complex Sq and an enzyme 5*2. Each reversible reaction is written as two 
irreversible reactions, e.g. the reactions ri and r2 are considered two distinct reactions and 
not one reversible reaction. 

2.3. Mass-action kinetics. Let J\f = {S,C,1Z) be a network. We denote the molar con- 
centration of species Si at time t by Cj = Ci{t) and associate with any complex y = 
{yi, . . . , yn) the monomial = HILi ^ For example, if y = (2, 1, 0, 1) G 

N , then the associated monomial is = cfc2C4. 

A rate vector is an element k = {ky^yi)y^yi G given by the assignment of a positive 
rate constant ky^y' G M+ to each reaction r: y ^ y' £ TZ. The (mass-action) species 
formation rate function corresponding to the rate vector k is defined by: 

M" ^ M", Yl ky^y'C^y' -y)- 

Let fK,i{c) denote the i-th entry of /«;(c), that is /K,i(c) = J^y^y'^n ^y^y'^^iv'i ~ Vi)- 

The set of reactions together with a rate vector give rise to a polynomial system of 
ordinary differential equations (ODEs): 

c = /.(c). 

These ODEs describe the dynamics of the concentrations Ci in time. The steady states of 
the network with rate vector k are the solutions to the system of polynomial equations in 
ci, . . . ,Cn obtained by setting the derivatives of the concentrations to zero: 

=/k(c). 

This system of equations is referred to as the steady-state equations. We are interested in 
the positive steady states, that is, the solutions c to the steady-state equations such that 
all concentrations are positive, c G W}. 
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Example 2.3. The ODEs system of the futile cycle taken with mass- action kinetics is: 

ci = -/C1C1C3 -I- {k2 + k3)c5 C4 = -^40204 + /cscs + fcece 

C2 = -k4C2C4 + (/C5 + kG)c6 C5 = kidCs - {k2 + kz)c^ 

C3 = -fcicics + A;2C5 + /cece c'e = ^40204 - (^5 + /c6)c6 

where the rate constant of reaction rj is denoted by ki. 

Remark 2.4. If k G and/or c G M" are not fixed then the function /k(c) can be seen 
as a polynomial function taking values in M[c], M[k] or M[cU k]. 

3. StOICHIOMETRICALLY COMPATIBLE STEADY STATES 

The dynamics of a network might preserve quantities that remain constant over time. 
If this is the case, the dynamics takes place in a proper invariant subspace of M", fixed by 
the initial concentrations Ci(0) of the system. Let v-v' denote the Euclidian scalar product 
of two vectors v,v' . Let f* denote the transpose of a vector v. 

Definition 3.1. The stoichiometric subspace of a network J\f = [S,C,TZ) is the following 

vector subspace of M": 

r = {y'-y\y^y'en). 

Two vectors c, c' G M" are called stoichiometrically compatible \i c — c' G F, or equivalently, 
if a; • c = a; • c' for all cj G F"*" . 

We denote by s the dimension of F. Note that F is independent of the choice of rate 
constants and depends only on the structure of the reactions. Being stoichiometrically 
compatible is an equivalence relation which partitions M_|^ into classes, called stoichiometric 
classes. In particular, the stoichiometric class of a concentration vector c G M" is {c + 

F}nl!^. 

For any rate vector k G MJ, the image of is contained in F. Thus, for any choice 
of rate vector k, c G M", and w G F-*-, we have that w • /k(c) = and thus a; • c = 0. 
If u; = (Ai, . . . , A„) then Ylll=i ^i^i independent of time and determined by the initial 
concentrations of the system. These conserved quantities are generally referred to as total 
amounts. Since a — 6 G F if and only if a; • a = a; • 6 for all oj G F-*-, total amounts are 
conserved within each stoichiometric class and characterize the class. 



Example 3.2. The stoichiometric space of the futile cycle (Example 2.2) is 

F = {Si + 53 — S's, S'5 — Si — 5*4, S'2 + 5*4 — 5*6, 5*6 — 5*2 — 5*3) 
= {Si + 5*3 — 55, 55 — Si — 54, 52 + S'4 — 5*6) 

and has dimension s = 3. The dimension of the orthogonal space F"*" isd = n — s = 3as 
well and a basis is F"*" = {uj^ ^uJ^ ,uj'^) with 

(3.3) a;i = 5i + S5, = 52 + Sg, = 53 + 54 + 55 + 56. 

Indeed, we have ci + C5 = C2 + ce = C3 + C4 + C5 + ce = 0. In this example, three total 
amounts determine each stoichiometric class. 
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Remark 3.4. Questions like "How many steady states does a system possess?" refer to the 
number of steady states within each stoichiometric class. If this restriction is not imposed 
and s = dim(r) < n, then the steady states describe an algebraic variety of dimension at 
least one over the complex numbers. 

Definition 3.5. A network M = (S,C, TZ) has the capacity for multiple positive steady 
states if there exists a rate vector k G and distinct vectors a, 6 G R" such that a — 6 E F 
and f^{a) = ff,{b) = 0. 

Note that in this work we focus mainly on the existence of multiple positive steady 
states. However, as we will show, the methods developed can preclude the existence of a 
specific type of multiple steady states on the boundary of . 

If a basis {w^, . . . , uj'^} of F-*- is chosen, then a network M has the capacity for multiple 
positive steady states if there exists a rate vector k £ MJ^ and distinct a,b £ M" such 
that /^(a) = = and cj* • a = cj* • 6 for all z = 1, . . . , In particular, if the map 

7^ : M" ^ Rd+n defined by /^(c) = {u^ ■ c, . . . ■ c, A,i(c), . . . , /«,„(c)) is injective, then 
the network does not have the capacity for multiple positive steady states. This function 
is the focus of study in this paper. 

If species in- or outflow reactions exist for all species in a network then the stoichiometric 
space has maximal dimension n. Therefore, the requirement a — 6 S F in the previous 
definition is superfluous. Specifically, if there is a reaction 5" — )• or — )• 5 for some 
species S then there are no vectors in F-*- with non-zero i-th entry. 

Lemma 3.6. Let M = {S,C,TZ) he a network. If Si ^ £ TZ or ^ Si £ TZ for some 
Si G S, then Aj = for all uj = (Ai, . . . , A„) G F-*-. 

Proof. Since ±5j G F, we have = w • (±5i) = ±Aj. □ 

Definition 3.7. Let M = {S,C,Tl) be a network. We say that 

(i) J\f is open if F = M". 

(ii) M is fully open if the outflow reaction Si ^ belongs to TZ for all Si £ S. 

(iii) M is closed if F 7^ M". 

If J\f is any network then the associated fully open network, J\f° = {S,C", TZ°), is the network 
with C" = C U S U {0} and n° = TZU {Si ^ 0\i = I, ... , n}. 

Lemma |3.6| ensures that a fully open network is also open. Fully open networks are 
considered by Craciun and Feinberg in [3] in the context of continuous flow stirred tank 
reactors and their results extend to arbitrary open networks. A closed network is allowed 
to have outflow reactions as well, but not for all species since F 7^ R". 

Remark 3.8. A closed network could also be made open by adding species inflow reactions 
to the set of reactions, but to be fully open requires species outflow reactions. 

Notation. Let M = {S,C,TZ) be a network and k a rate vector. At this point we have 
defined the species formation rate function J^, the stoichiometric space F, and used n for 
the number of species, s for the dimension of F and d = n — s for the dimension of F-*-. 
None of these objects incorporate reference to the specific network in the notation. This 
notation is used without further mentioning throughout the paper. Additionally, to ease 
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notation, in some examples species are called A, B, . . . and the respective concentrations 

CA,CB, 

4. Degenerate steady states 

For any function / = (/i, . . . , fp) : — )• W let Jc(/) denote the Jacobian of / at c, 
that is, the p x m matrix with entry (i, j) being dfi{c)/dcj. 

Definition 4.1. Let N = (5, C, IZ) be a network and k G a rate vector. A steady state 
c* G M" of is degenerate if ker(Jc*(/«)) n T / {0}. 

It is proven in [9, §5] that for any n G M^, c G M" , and 7 G M", 

n 

(4.2) Jc(/«)(7) = ^ ky^y,cy{y *cl){y' - y), where u *c w; = ^ ^-^ . 

y^y'en i=l * 

It is our aim to understand ker( Jc*(/k)) H F. We find a criterion to determine whether 
this intersection consists of the zero vector or is a proper subspace. The criterion does 
not require the computation of ker(Jc* (/«)). In particular, we find that the existence of 
degenerate steady states is linked to the function being injective. 

Since vectors of F are characterized by being orthogonal to all vectors in F-*- , we obtain 
the following proposition (stated here for a general vector subspace F). 

Proposition 4.3. Let Af be a network and k G a rate vector. Let F be a vector 
subspace o/M", . . . ,0;™} a basis of F^ and define : M" R'^+" by 

1k{c) = (w^ • c, . . . ,0;'" • c,/k,i(c), . . . ,/k,„(c)). 

Fix c* G M". Then, ker( Jc.(/k)) n F = {0} if and only if the Jacobian Jc-if^) of f 1^ at c* 
has maximal rank n. 

Proof. Let Jc*{f^) be the Jacobian of /^^(c) at c* . It is an (m + n) x n matrix. The rank 
of Jc*{f k) is maximal if and only if ker( Jc*(/^)) = {0}. For i = 1, . . . ,m, the z-th row of 
Jc*(/^) equals the vector w*. The lower n x n matrix of Jc*{f n) (obtained by removing 
the first m rows) is equal to Jc*{fK)- A vector v belongs to F if and only if • u = for 
all i. It follows that v G ker( Jc*(/k)) if and only if w G ker( Jc*(/k)) and v £ F. Thus, 
ker( J,* (/,)) n F = {0} if and only if ker( J^. (/J) = {0}. □ 

By letting F = F in the previous lemma, we have shown that a steady state c* is non- 
degenerate if and only if Jc'ifn) has maximal rank n. Since Jc*(/k) is a (n + d) x n 
matrix, d rows of the matrix are linearly dependent of the remaining n. We describe now 
a procedure that is independent of c and k to determine d rows with this property. Thus, 
it can be used to determine the rank of Jc*{ff^). 

Lemma 4.4. If 00 £ T-^ then oJ^Jcifn) = for all c G M" and all rate vectors k G MJ. 

Proof. If cj G F"*" then a; • /^(c) =0 for all c and hence the scalar product vanishes as 
a polynomial in c. It follows that = '^^'^^"^^^^ = <^ " ^^g}^^ for all i = 1 , . . . , re and 

Uj'Jc{fn)=0. □ 

Hence, each io G F-*- provides a vanishing linear combination of the rows of Jcifn)- 
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Definition 4.5. A basis {u!^, ■ ■ ■ , uj'^} of F-*" with = (A^, . . . , A^) is said to be reduced if 
A* = 1 for all i and A*- = for all j = 1, . . . , j, . . . , d. Given a rate vector k, the associated 

extended rate function : M" — t- M" is the function defined by 

/k(c) = {iO^ ■ C, . . . ,Uj'^ ■ C, f^^d+l{c), . . ■ , fK,n{<^))- 

After reordering the species, such a basis of F"*" always exist (use for instance Gaussian 
elimination on any basis of r-*"). We assume from now on that a reduced basis is chosen, 
implying that the species are ordered accordingly. 



Example 4.6. The basis uj^,uj'^,uj'^ of of the futile cycle given in (3.3) is reduced. The 
associated extended rate function is 

= (Cl + C5, C2 + C6, C3 + C4 + C5 + Cg, -/C4C2C4 + ksC^ + k^CQ, 

klCiC3 - {k2 + k3)c5, A;4C2C4 - (fcs + A;6)C6)- 

Let Jciifn) denote the i-th. row of Jc Hk)- If W^, ■ ■ ■ > w'^} is a reduced basis with = 
(A^, . . . , A^), then it follows from Lemma |4.4| that 







(a;*)*Je(A) = Jc,.(/«) + Yl ^PcAh)^ 

j=d+i 

Thus, the rows 1, . . . , d of Jdfn) are linear combinations of the rows d + 1, . . . ,n. It follows 
that the rank of Jc*^/^) equals the rank of Jc*(/k)- The latter is a square n x n matrix 
and has maximal rank if and only if its determinant does not vanish. 

Corollary 4.7. Let Af be a network, k E a rate vector, {co^, . . . jUj"^} a reduced basis of 
r-*- and /k(c) the associated extended rate function. If c* G M", then ker( Jc* (/«)) nL = {0} 
if and only i/ det( Jc* (/«)) 7^ 0. In particular, a steady state c* € M** of J\f is degenerate if 
and only i/ det( Jc* (/«)) = 0. 

Example 4.8. The Jacobian matrix Jcifn) of the futile cycle associated with the reduced 
basis of in (3.3) is 



Ufn) 



The determinant of Jcifn) is 

(4.9) - det(Jc(/K)) = (ciC2 + ciC4)kik3k4, + cikik^k^ 

+ (ciC2 + C2C3)kik/^k^ + C2k2kik^ + C2k3kik^ + cikik^k^. 

All coefficients of —det( Jc(/k)) as a polynomial in c,k are positive. Thus det( Jc(/k)) 
does not vanish for any c G M" and n G and all positive steady states of the futile 
cycle are non-degenerate. In this example, det(Jc(/K)) is linear in the rate constants. If 
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at least one coefficient in tlie polynomial was negative then we could find c, k for which 
det(Jc(/K)) = 0, implying that a degenerate steady state could occur. This observation 
holds for any network and will be discussed in the following sections. 

Remark 4.10. The minimal space containing the image of is the kinetic suh space: 
(4.11) A, := (im/,) cr. 

In general, the two spaces and T might not agree for a fixed rate vector k. If this is the 
case then ker( Jc(/k)) R F 7^ {0} for any c G M , : as above, if u; G Kj: then = a; • = Ci; • c. 



Similarly to the proof of Lemma 4.4 we have that io^Jcifn) = and hence (im Jcifn)) ^ ^k- 
If C r then {im Jcifn)) has at most dimension s — 1 and it follows that the dimension 
of ker (Jc(/k)) H F is at least 1. Thus, if the stoichiometric and the kinetic spaces do not 
agree for some k, then all steady states corresponding to are degenerate. 

Example 4.12. Consider the network with reactions ri: A ^ B and r2 : A — )• C. The 
species formation rate function : — )• M'^ is 

/k(ca, cb, cc) = {-{ki + k2)cA, kiCA, k2CA), 

where ki denotes the rate constant of reaction rj. The stoichiometric space has dimension 
2 and a basis of F-*- is A + B + C. However, k2B — kiC G A;^ \ F"*- and thus for all k, all 
steady states are degenerate. 

Example 4.13. The stoichiometric and the kinetic spaces of the previous example never 
agree. For some networks, the two spaces do not agree only for certain rate vectors. 
Consider for example the network with reactions ri : A ^ B , r2 '■ A ^ C and rs : B + C 
2 A. The species formation rate function : — )• is 

/k(ca, cb,cc) = {-{ki + k2)cA + 2k3CBCc, kiCA - ksCBCc, k2CA - ksCBCc). 

If we let uj = (Ai, A2, A3) then the equation oj ■ fnic) = for all c gives 2Ai = A2 + A3 and 
(A2 — A3)(A;i — k2) = 0. If ki / k2, then (1, 1, 1) is a basis of F-*- as well as A;^ and the 
stoichiometric and the kinetic spaces agree. However, if ki = k2 then (1,2,0) G A;^ \ F-*- 
and the two spaces do not agree. Further, in this case all steady states are degenerate. 



Remark 4.14. Note that Proposition 4.3, Lemma 4.4 and Corollary |4.7|do not depend on 



the kinetics being of mass-action type, but hold for general (differentiable) kinetics with 



and /k changed accordingly. In particular. Corollary 4.7 characterizes degenerate steady 



states in terms of the determinant of the Jacobian of /«. 

5. InJECTIVITY OF CHEMICAL REACTION NETWORKS 

Here we introduce the notion of injectivity: a network is injective if for all rate vectors 
K the function is injective over M" . The definition is an extension of the definition of 
injectivity for fully open networks j4i Def. 2.8] to arbitrary networks and, together with 
Proposition |5.2| below, it is in line with the definition given in a recent paper on concordant 
networks [17J. We show that a network is injective if and only if the Jacobian of is 
non-singular for all positive concentrations c and for all k . 

In [Tjj Def. 6.1] a definition of injectivity is given for networks with F C M". It relies 
on the Jacobian of the species formation rate function restricted to the stoichiometric 
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space. We prove in Theorem 5.6 below that our definition agrees with their definition. The 
equivalence is also claimed in Remark 6.4 in |^ without a proof. 

Definition 5.1. A network M is said to be injective if for any distinct a,b G M" such that 
a - 6 G r, we have /^(a) / 

The definition of injectivity is restricted to the interior of the positive orthant. However, 
we show below that being injective also precludes the existence of distinct a, b such that 
a — b G T and /^(a) = fnib), provided either that a G M_,_, b G M" or that a, 6 G are 
on the boundary of and fulfill a certain condition. For a G , let la = {i\a.i = 0} 
be the indices for which Oj is zero and let 3^a = {y\y ^ u' ^ T^i la H supp(y) ^ 0}, where 
supp(y) = {i\yi / 0} is the support of y. That is, ya is the set of reactant complexes 
involving at least one species Si for which Oj = 0. It follows that y belongs to 3^a if and 
only if = 0. 

Proposition 5.2. Let M be a network. The following two statements are equivalent: 
(i) N is injective. 

(a) For any distinct a,b £ such that 3^a H 3^;, = and a — b GT, /^(a) 7^ fnib). 
In particular, if b G M" , that is = 0, then 3^a H 3^;, = is fulfilled for any a G M^. 

Proof, (ii) obviously implies (i). To prove the reverse we assume that there are a,b G M", 
such that 3^a n 3^;, = 0, 6 — a G r and /^(a) = /k(^) and show that A/" cannot be injective. 
Let 7j = 6i — Oi and for a set M C {1, . . . , n} and y G C define x^^^ = HieAf Further, 
define I = laD lb and J = I^. If / = 0, then clearly J\f cannot be injective. Hence, assume 
that / 7^ 0. We seek to define a rate vector k = {ky^y')y^y' and a,b such that di = Oj, 
bi = bi for i /, di, bi > and 6j — Oj = 7^ for i G /, and 

~ky^y.{dy -by){y' -y)=Q, 

that is, fhid) = fh{b). Then, since d, 6 G M" and 6 — a G F, it follows that M is not 
injective. The equality /^(a) — = can be rewritten as: 

0= E ky^y'{ay^by){y'-y)+ ^ W " " y)> 

y^y' y^y' y^^y 

yey^^ny;; yey^nyi; yey^^ny^ 

Assume that we can find dj, 6j > with bi — di = 7^ for all i, di = Oi, bi = bi for i ^ I, and 
ay > such that 

(5.3) — y--by = dy, iiyGyany^, and — ^as' = 6^ if y G 3^,^ n j;^. 

ay + 1 ay + 1 

Then, let fc^^y/ = fcj^^j^/ for y G 3^^ n y§, ky^y, = {ay + l)ky^y/by/by for y G 3^a n 3^^ and 
/cy_!.y' = {ay + l)ky^yiay /dy foic y G n J'ft. With these definitions, if y G 3^a H 3^^ since 
= we have 

ky^y'{dy - by) = ky^y,{ay + l)(d^&V6^ " 6^) = ky^y,{ayby " (Oj, + 1 ) 6^ ) 
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Analogously, ky^yi{ay — by) = ky^yi{ay — h^) if y G 3^^ n y^,, using = 0. It follows that 
for k = {ky^y/)y^y/, WB havc fii{d) - fkib) = f^{a) - f^,{b) = as desired. 



All that remains is to find di,bi > for i £ I and > such that (5.3) is fulfilled. If 



y ^ 3^6) then by,dy are independent of the value of di,bi for i G If,. Therefore, we define 



first di,bi for i ^ la and focus on fulfilling (5.3) for y G 3^a H yp. Further, since we want 



bi — di = 7j, once a, is defined, we have 6j = Oj + 7j. 

Let /a = {?!,... , im}, Ik = {Ht ■ ■ , ik} {k < m) with Jq = and Jk = J U Ik- Further, 
recursively define = {y G 3^a| la H supp(y) C \ y^-i with yo = 0. It follows that 
UfcLi = 3^ such that any complex in y is in precisely one 3^^. That is, if y G J'fc then 
ik is the largest index in in the support of y. We define di^,bi^ recursively. Note that 
7i > for all i G /«• Assume that Oj^, . . . , aij._i and bi-^, . . . , bi^_^ are defined and hence 
the products ft^-^^-i, a^-^^-i are fixed. For Zfc do the following: if yk = 0, let di^ = 1 and 



bi^ = l + 7jj.. If 3^fc 7^ 0, observe that for y G yk (for which yi^, ^ 0), the equality (5.3) can 
we written as 

^ + l^kr'' " % + 1 a^^-^ ' l^-r- (d- V ) ' 

with ri^{di^^,yi^) := df^''/{di^ + HkY''' > 0. The function r^^ is increasing in ai^ and 
Tj^, (0,yj^,) = for all y^^. By defining di^ > arbitrarily such that 

' lyJk-i 



X^^ik^yik) < mm : 



'^ik v^jfc ! tf «fc ; ^ ^^^^^^ \^ 



y G 3^fc 



for all y ^yk (which is a finite set), we obtain > fulfilling (5.4) as desired. 

The same procedure is applied to define Oj, bi for i £ h, with the roles of a and b reversed. 
In this case, however, 6j = implies that 7^ < and di = bi — ji for i G It- Therefore, 



^iki^ikiVik) becomes b-'^'' /{bi^ - 74^)^'*, which also is increasing. □ 

The assumption J'a H 3^;, = in Proposition |5.2[ ii) cannot be relaxed. Consider the 
network with 5 = {A, B},C = {A + B,A + 2B, 0} and reactions A + B ^ 0, A + 2B ^ 0. 
This network is injective and open. If a = (1, 0) and b = (0, 1) then 3^a H 3^^ = {A + B, A + 
2B} 7^ and for any rate vector k we have /^(a) = fnib). 

Remark 5.5. In [17], a network M with arbitrary kinetics is said to be injective if /«(«) 7^ 
/k(6) for any pair of stoichiometrically compatible concentration vectors 0,6, at least one 



of which is positive. The condition given in Proposition 5.2 ^ii) is slightly more general in 
that both a and b can be non-negative. 

If M is injective, then for any choice of rate vector k at most one positive steady state 
can exist within each stoichiometric class, i.e., for every cq G M" there exists at most one 
c G {co + r} n such that /k(c) = 0. In other words, if J\f is injective then J\f does 
not have the capacity for multiple positive steady states. However, the reverse might not 
be true: non-injective networks exist that do not have the capacity for multiple positive 



steady states. An example is provided in Example 7.8 



The proof of the following theorem is adapted from the proof of [4^ Th 3.1]. 
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Theorem 5.6. Let M = {S,C,Tl) be a network. Then, M is injective if and only if 
ker( Jc(/«)) n r = {0} for all cGRI and k £ 



Proof. N is not injective if and only if there exists n G and distinct a,b £ M" such that 
a - 6 G r and f^{a) = f^,{h). Further, ker(Jc(/^)) n T / {0} for some c G M*;: and r/ G 
if and only if there exists 7 G T such that Jc{fri){l) = 0. By definition and using (4.2), 

We will show that given distinct a,b €z M" such that 7 := a — 6 G T and some rate vector 
K, there exist c G M" and rj G such that 

(5.7) ky^y'iay - by) = r]y^y'cy{y *c 7) 

and vice versa. Consider distinct a,b €z M" such that 7 := a — 6 G T, 7 7^ 0. If Oj — 6i 7^ 
define Cj := > and let Cj = 1 otherwise. Since 7^ = if = 6j, we have 

ay = byey*"^ . Note that the signs of ay — by = by{ey*''^ — 1) and y*cl agree. liy*cl = 0, let 
rjy^yi = 1. Otherwise, we let r]y^yi = ^~^y(^y^ — > and (5.7) is fulfilled. Reciprocally, 
given 7 G r, 7 / 0, c G M!j: and r/ G MJ, define a, 6 G M!J: by 6, = = 1 if 7^ = 
and bi = ^i/{e'^^^^^ — 1) > 0, Oj = bie"'^/^^ > otherwise. Then a — 6 = 7 G F and 
gV — }jy = i,y(^gy*ci _ j^g^g ^j^g same sign as y *c 7. If ay — by = 0, define ky^yi = 1. 

Otherwise, define ky^yi = ^l^i^av-bv)^ ^ ^ ^^"^ equality (5.7) is fulfilled. 

□ 

Remark 5.8. The summand in corresponding to an infiow reaction — )• y' in TZ takes 
the form kQ_^yiy' and thus is independent of the concentration vector c. It follows that 
<^c(/k) is independent of the presence or absence of infiow reactions in TZ and so is the 
property of being injective. 

We have thus obtained a characterization of injective networks in terms of the Jacobian 
associated with the species formation rate function. Together with Corollary 4.7 we obtain: 

Corollary 5.9. Let Af be a network, {oj^ , . . . , uj'^} a reduced basis of and /k(c) the 
associated extended rate function. N is injective if and only if det(Jc(/K)) 7^ for all 
K G M5 andc£ Wl. 



Remark 5.10. Statements similar to Theorem 5.6 and Corollary 5.9 cannot be stated for 
individual k, since k and 77 generally are different in the proof of Theorem 5.6 However, if 
the total degree of each of polynomials in the components in is at most two then is 
an injective function if and only if the Jacobian is non-singular [3]. 

Note that M is injective if and only if the extended rate function associated with 
a reduced basis is injective over M". Further, as is observed in (4.9) for the futile cycle, 
det(Jc(/K)) is a homogeneous polynomial in the entries of k with total degree s, because 
the rows 1, . . . , d of Jc(/k) are constants (that is, independent of c and k). 
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Remark 5.11. A general version of the corollary above has recently been formulated in 
[13j for a certain class of polynomial maps, without the restriction to species formation 
rate functions and conservation laws associated to networks. 



From Definition 4.1 and Theorem 5.6 we obtain the following corollary. 



Corollary 5.12. Let M be a network. If M is injective then there exist no degenerate 
positive steady states. 



Remark 5.13. It follows from Theorem 5.6 and Remark 4.10 that if the two spaces F and 



Ak are not identical for some rate vector k then the network cannot be injective. 

As noticed in Example |4.9[ the determinant of the Jacobian of the extended species rate 
formation function of the futile cycle can never vanish. Thus, the futile cycle is injective 
and does not have the capacity for multiple positive steady states. 

The coefficients of the determinant of the Jacobian of fully open networks are charac- 
terized by Craciun and Feinberg in [4J and this characterization easily generalizes to open 
networks. Thus, in order to characterize the coefficients of the determinant of the Jaco- 
bian of an arbitrary network M, we consider the associated fully open network, J\f°, and 
"match" the terms of the respective determinants. This is done in Rafter we discuss some 
results about open networks in the next section. 

6. INJECTIVE OPEN CHEMICAL REACTION NETWORKS 
Recall that a networ k is open if its stoichiometric space is M". If this is the case then 



/k = /k- By Theorem 5.6 an open network M is injective if and only if the Jacobian 



Jcifn) is non-singular, i.e., det(Jc(/K)) 7^ for all rate vectors k G and all c G 



Hidden in the proof of Theorem 3.1 in [4J and Theorem 5.6 above we find a simplification 
of the characterization of injective open networks: for a network to be injective it suffices 
to fix any concentration vector c G M" and show that <let{Jc{fn)) 7^ 0, for all rate vectors 
K G M^. 

Proposition 6.1. An open network N is injective if and only if det{Jc{fK)) 7^ for all 
rate vectors k G and a fixed positive c G WL ■ 



Proof. By Theorem 5.6, it is enough to prove that det( Jc(/k)) 7^ for all rate vectors 
K G and all c G if and only if the statement holds for a fixed c G . The 
forward implication is obvious. To see the reverse, assume that det(Ja;(/K)) = for some 
rate vector k G and x G M", that is, there exists a non-zero vector 7 G such 
that = Jx(/re)(7)- Define rj = {r]y^yi)y^yi G with r]y^yi = ky^yixy /c^ > 0, and 
5 = . . . , 5n) with 5i = '^iCi/xi. With these definitions, we have ky^y/x--' = rjy^yic^ and 
y *x 1 = y *c ^- Then, using (|4.2|), we have 



Jc{fv){^) = E *c S){y' -y) = Ufn){l) = 0, 

which implies that det(Jc(/,j)) = and we have reached a contradiction. □ 

The simplification presented here applies to any open network, independently of whether 
it contains outfiow reactions or not. We have proved that injectivity of an open network 
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can be checked using the Jacobian criterion with the concentration vector fixed to 1 := 
(1, . . . , 1). In this case the determinant of the Jacobian is a polynomial depending only 
on the rate constants and this reduces the number of variables substantially. Further, the 
polynomial is homogeneous of total degree n. 

In [H Rk. 2.9] a different simplification is performed where the rate constants of all 
species outflow reactions are fixed to 1. We state this result as a proposition below and 
give the proof for completeness. 

Proposition 6.2. LetM be an open network and fix a subset O C {1, . . . , n}. The following 
statements are equivalent: 

(i) M is injective. 

(ii) det(Jc(/K)) 7^ for all c G M" and all rate vectors k satisfying ks^^o = 1 for all 
i£0. 

(Hi) For any distinct a,b £ M", 7^ /k(^) for all rate vectors k satisfying ks^^o = 1 

for all i G O. 



Proof, (i) implies (ii) and (iii) by Definition 5.1 and Theorem 5.6 We now prove that (ii 



and (iii) separately imply (i). For any k G and z G M" define k = {ky^yi)y^yi G 
and z G M" by: 

^_\ziks^^Q, i£0 Z _ ^y^y' 

K,,, ' 



Zi I 



Clearly ks-^o = 1 for all i G O, Ky^y'(a)y = ky^y/a^, and further y 7 = y 7 for 
any 7 G M". It follows that /^(a) = /«(«), and, similarly, for any 7 G M", Juifn)^) = 
Ja{fK){l)- Therefore, if M is not injective, then (1) there exists a rate vector k G and 
distinct a, 6 G such that /^(a) = /k(5) and the construction of k G and o, 6 G 
as above implies that (iii) does not hold, and (2) there exist a rate vector k G M^, a G W\. 
and 7 G M", 7 / 0, such that Ja{fK){l) = and thus a rate vector k G M^, a G 
and 7 G M" contradicting (ii). It follows that (ii) and (iii) both imply (i) and the proof is 
completed. □ 



The simplifications in Propositions 6.1 and 6.2 cannot be performed at the same time 
because it would constrain the choice of the free variables k and c too much. 

The next proposition is an extension of [H Eqn. 3.15] where the statement is made for 
a specific class of open networks. However, the proof works line by line for the class of 
all open networks. In fact it does not depend on whether the network is open or not, but 
all terms in the proposition are zero unless the network is open. Notice however that our 
statement differs from the statement by Craciun and Feinberg in [3] in the sign (— 1)", 
because they establish the Jacobian criterion on the Jacobian of — /«. 

Recall that n is the number of species and note that det(Jc(/K)) is a homogeneous 
polynomial of degree n in the rate constants. To simplify the notation we introduce the 
following: for any set of m reactions, R = {y^ y'^ , . . . , y"^ — )• y""}, let 

• y{R) be the n x m matrix whose i-th column is y*. 

• r(i?) be the n x m matrix whose i-th column is y* — y'*. 

• If m = n, define a{R) = (-1)" det(3^(i?)) det(r(i?)). 
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Proposition 6.3 ([3]). Let M be an open network and let R = {y^ — )• y'^, . . . , — )• y'"} 
be a set of n reactions. Viewed as a polynomial in the rate vector k, the coefficient of the 
monomial YYl^i kyi^yii in (iet{J c{f k)) with c G M" is 

a{R) :=a(i?)c-^+^-i^\ 
In particular, the coefficient of the monomial Y\!i=i kyi^yH in det(Ji(/K)) is cr{R). 

Remark 6.4. The term a{R) is a monomial in c: if Y17=i u] ~ 0' then Sj has zero 
coefficient in for aU i and thus = det(y^, . . . , y") = a{R). If a reaction — )• y" 
appears twice in a set R then cr{R) = 0. Therefore the degree of each ky^yi in the 
polynomial det(Jc(/K)) is either zero or one. 

Proposition 6.5. An open network M is injective if and only if the non-zero coefficients 
(T(i?) have the same sign for all sets R of n reactions and there exists at least one set for 
which (j{R) 7^ 0. 



Proof. The reverse implication follows from Corollary 5.9, For the forward implication, 
assume that M is injective. Clearly, det(Ji(/K)) / and thus there exists at least one set 
R for which (y{R) 7^ 0. Note that det(Ji(/K)) is a polynomial in k with total degree n 
and degree at most one in each variable ky^yi. Assume that there exist two coefficients 
Ri,R2 satisfying a{Ri) > and a{R2) < 0. Set ky^y' = 1 if y — )• y' ^ i?i and ky^y/ = k 
if y — )• y' G where is a positive parameter. After this transformation, the monomials 
corresponding to sets of reactions R ^ Ri have degree in k strictly lower than n. Then 
det( Ji(/k)) is a polynomial of degree n in A;, with positive leading coefficient. It follows that 
if k tends to +00, then det(Ji(/K)) > 0. Symmetrically, using R2 we can find rate constants 
for which det(Ji(/K)) < 0. Since det(Ji(/K)) is continuous in a connected domain, there 



exists a rate vector for which det(Ji(/K)) = 0, contradicting Proposition 6.1 □ 



The criterion in Proposition |6.5| is independent of the rate vector k. The requirement 
that there exists at least one set R for which cr{R) 7^ cannot be removed. Consider for 
example the network M with set of reactions {A — )■ B,A^ 0}. The stoichiometric space 
has dimension 2 and thus M is open. However, det( Jc(/k)) = for all k,c, since cb is 
not a variable of This requirement is not imposed in jU Th. 3.3] because it holds 
automatically for fully open networks: in fact, the set of reactions {Si — >• 0, . . . , 5„ — >• 0} 
provides the non-zero coefficient (— 1)" in the determinant expansion of the Jacobian. It is 
mentioned in [H Rk 3.5] that the requirement is necessary if the network is not fully open. 

7. Injective closed chemical reaction networks 
We would like to have a characterization of the coefficients of the polynomial det( Jc(/k)) 



for closed networks similar to that of Proposition 6.3 and a characterization of injectivity 



similar to that of Proposition 6.5 , To this end we consider the fully open network associated 
to a network and use the results of the previous section. 

Before getting into technicalities we illustrate the idea with the futile cycle. The futile 
cycle has no outflow reactions and therefore the fully open network has an extra reaction 
Si ^ for all i. Let g^o denote the species formation rate function of the fully open futile 
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cycle with any rate vector such that ks^^o = 1- Then, the system of ODEs of the fully 
open futile cycle taken with mass-action kinetics is: 

C'l = -kiCiCs + (/C2 + k3)c5 - Ci 

c's = -kicics + k2C5 + keCQ - 03 



C"g = kiCiC3 - {k2 + k3)c5 - C5 
The Jacobian Jc{gK°) of 5^° is 

/ -k}C3 - 1 



C2 

C4 

C6 



-/C4C2C4 + (fcs + /C6)C6 - C2 

— /C4C2C4 + A;3C5 + A;6C6 — C4 
/C4C2C4 - (A;5 + ke)cQ - cq. 



V 



-fclC3 



-kiC3 


A;iC3 






-A;4C4 — 


—^404 


A;4C4 



1 



-kici 


kiCi — 1 








— /C4C2 


-k4^C2 - 1 



fc4C2 



A;2 + fcs 

A:2 

-h - ks- 






fes + ^6 

fee 

^6 



ke - 



\ 



1 



Propositions 6.2 and 6.3 link the properties of the coefficients of the determinant of Jc{gK°) 
to the injectivity of the fully open futile cycle. The determinant of Jc{gK°) does not change 
if the fifth row is added to the first, the sixth row to the second, and the fourth, fifth and 
sixth to the third. Thus, 



det(Jc(5'Ko)) 



where the sign — in front corresponds to changing the sign of the first three rows. This 
determinant is almost equal to the determinant of Jcifn) as one can see from Example 



1 











1 








1 











1 








1 


1 


1 


1 





— k^Ci 





-A;4C2 - 1 


^3 


h 


kiC3 





felCl 





-k2 - ks-l 








A:4C4 





A;4C2 





-fcs - ke 



4.8 



The difference between the two determinants arises from the —1 in the diagonal entries 
of the matrix for the rows 4,5,6. Therefore, by splitting column 4 using (0, 0, 1, — A;4C2 — 
l,0,/c4C2) = (0, 0, 1, — A;4C2, 0, A;4C2) + (0,0,0,-1,0,0), and similarly for columns 5,6, we 
have 

det( Jc((7reo)) = — det( Jc(/k)) + monomials of total degree at most 2 in ki, kQ. 

The determinant det( Jc(/k)) is a homogeneous polynomial in ki, . . . ,kQ of degree s = 3. 
Thus, it agrees with the terms in — det{Jc{gK°)) of total degree 3: 

- det{JcigK.°)) = (ciC2 + ciCA)kik3ki + ciikiksk^ + kik3k^) + (ciC2 + C2C3)kikik^ 

+ C2{k2k4,k^ + ksk^k^) + (ci + C3)(A;iA;6 + kik^) + ^2^6 + ^3^6 + ^2^5 

+ ^3^5 + C2k/^k^ + (C2 + C4,){k2k4 + kski) + (ci + C3)(C2 + C4)kik4 

+ cikiks + (ci + C3)ki + k2 + k3 + (c2 + 04)^4 + k^ + k^ + 1. 

The fully open network is injective since all coefficients of minus the determinant expansion 
are positive. 

In the example, the row modifications done prior to the computation of the determinant 
were based on the fact that the futile cycle has conservation laws. This principle holds 
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generally for any closed network. Let N" be the fully open network associated with a 
closed network J\f = {S,C,Tl) with stoichiometric space F. Let 0{J\f) := {i\Si — )• ^ IZ] 
be the set of indices for which species outflow reactions do not belong to TZ. Then 1Z° = 
n\J{Si ^ 0| i G 0{N)}. For example, if TV is the futile cycle then 0{N) = {1, 2, 3, 4, 5, 6}. 



If {oj^ , . ■ ■ ,1^^^} is a reduced basis of F"*" then by Lemma 3.6 we have that i G 0{Af) for 
i = 1, . . . , d. Thus, the cardinality of 0{M) is at least d. 

We start by relating the species formation rate functions of M and J\f". Given a rate 
vector K = {ky^yi)y^yi S M J , define the associated rate vector k° = {k° i)y^yi G by 
setting ^5.^0 = 1 if « € 0{N) and k°^y, = ky^y/ for y ^ y' gU. 

Let 5i = 1 i £ 0{N) and 5j = otherwise, and let EJf^ be the n x n matrix with 
zeroes everywhere but 5i in the diagonal entries for i = m, . . . ,n. Then, the species 
formation rate function g^o (c) of M° with rate vector k° is 

5k°(c) = /k(c) - ((5lCl, . . . , 6nCn)- 

Theorem 7.1. Let M be a closed network, {ui^ , . . . , lo'^} a reduced basis ofT^ and fnic) 
the associated extended rate function. For a rate vector n £ M^, let n° , g^o^c), and 5i 
be defined as above. Then, 

detiMg.o)) = {-!)'' det{JciL)-E2+,). 

Proof. Since 5ko(c) = /^(c) - ((5iCi, . . .,5nCn), we have that Jcign") = Jcifn) - Ef. We let 
uj'^ = {X\, . . . , A^). If 5j = 0, then 5j — t- € 7^ and it follows from Lemma 3.6 that A* = 
for all i. Thus, {uj^YEf = {6iX\, . . . ,(5„A^) = and from Lemma 4.4 

(u'fMg.o) = {uJUfn) - (t^')*^r = for i = 1, . . . , d. 

Let P be the n x n matrix whose i-th. row is — for i = 1, . . . ,d and the i-th. canonical 
vector ej for i = d + 1, . . . , n. By the choice of w*, det(P) = (—1)'^. Further, PJc{gK°) = 
Jc(/k) - E'^+i and hence, 

{-if detiUg^o)) = det{PJc{g^o)) = det{ML) - E^+i). 

□ 

Corollary 7.2. The determinant expansion of Jcifn) o-s a polynomial in k agrees with the 
terms in the determinant expansion of {—if det{Jc{gK°)) of total degree s. 

Proof. For any non-empty set / C {d + I, . . . ,n}, let Jc{fK) be the matrix whose i-th 
column equals that of Jc(/k) for i ^ / and is the vector -SiCi for i £ L Then, by the 
column multilinear expansion of the determinant, we have 

(-l)'^det(J,(5..^)) = det( J,(7«))+ Yl det(4(7,)). 

0^7C{d+l,...,n} 

If c is fixed then det(Jc(/K)) is a homogeneous polynomial in k of total degree s, while the 
terms det(J^(/K)) are polynomials in k of total degree strictly lower than s. □ 

Let be a network. Let TZg be the set of all sets of n reactions formed by the union of 
a set R containing s reactions in TZ and d = n — s outflow reactions S'j — )• with i G 0{J\f). 



Since the cardinality of 0{M) is at least d, TZs / 0. By Corollary 7.2 and Proposition 6.3 
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the coefficients of det(Jc(/K)) as a polynomial in k and c are, up to a sign, a{R') for 



R' G IZ,. If i 



J' 



1, . . . ,d are the indices for which outflow reactions Si 



belong 



to then det(3^(ii')) and det(r(i?')) are simply the minors of y{R') and T{R!) with the 
ij-th rows and columns removed. Equivalently, det(3^(i?')) and det(r(i?')) are the minors 
obtained by removing the ij-th. rows from y{R) and r(i?), respectively, for R = R' D TZ. 



The matrices y{R) and T{R) are introduced on page 14 



Let Od{M) be the set of subsets of 0{N) of cardinality d. If M is any n x s matrix and 
/ G Od{M), let M/ denote the s x s submatrix of M obtained by removing the j-th row for 



all j G /. The following corollary is a consequence of Proposition |6.3| and the discussion 
above. 

Corollary 7.3. Let N he a network and s he the dimension of the stoichiometric space. 
The terms in the expansion of the determinant det{Jc{fK.)) are monomials in k of total 

..,y'^ y'^} be a 
in det{Jc{fK.)) for 



degree s and linear in each rate constant. Further, let R = {y^ 
set of s reactions from M . The coefficient of the monomial k. 

cGRI is 



Yl detiyiR)i)detiTiR)j) 



Alternatively, the coefficient of 



=1 ^y* 



can he written as 



a. 



R'eHs,R=R'mz 



Observe that the vector ^l^i is simply the row sum of the matrix y{R). The corollary 
reduces to Proposition |6.3| if J\f is open. 

Remark 7.4. Only the rate constants of outflow reactions that are not in J\f, i.e. those in 
0{J\f), are set to one in the associated rate vector k°. Otherwise the determinant of Jcifn) 
would not be a homogenous polynomial of total degree s in k. 



Example 7.5. Consider the futile cycle and (4.9). The coefficient of kik^k^ in — det( Jc(/k)) 
is ciC2 + C1C4. It corresponds to the reactions R = {ri,r3,r4} and the matrices y{R) and 
r{R) are: 

\ 

1 



y{R) 



( 1 


1 






1 


0/ 



r(i?) 



-1 



V 







-1 



-1 
1 




\ 

1 

1 


1 / 



The only sets of indices / G 0^{J\f) for which the product det(3^(ii)7-) det(r(i?)/) is non- 
zero are / = {1,2,6} and / = {1,4,6}. These sets give the coefficient C1C2 + C1C4. 
Since the last row of y{R) is zero, the index 6 belongs to all index sets / for which 
d(ii{y{R)i) det(r(i?)/) / 0. 



Corollary 7.6. Let Af 

(i) M is infective. 



{S,C,Tl) he a network. The following are equivalent: 
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(ii) The non-zero coefficients <7{R') have the same sign for all sets R' G IZg, and cr{R') / 

for at least one set R' . 
(Hi) The non-zero products det(3^(i2)/) det(r(i?)7-) have the same sign for all sets R of s 

reactions in TZ and I £ OdiM), and further det(3^(i?)/) det(r(i?)/) 7^ for at least 

one set R and some I G OdiM). 



Proof. The equivalence between (ii) and (iii) is a consequence of Corollary 7.3 If (ii) holds 



then (i) is a consequence of Corollary 7.3 and Corollary 5.9 To show that (i) implies 



(ii), we use the same argument as in the proof of Proposition 6.5 Using Corollary 7.3 and 



Corollary 5.9 it suffices to show that for any set R! G TZg we can find a concentration vector 
c and a rate vector k for which the sign of det(Jc(/K)) agrees with the sign of {—iya{R'). 
Let /c be a positive parameter and let m be an integer. Let R = R' CiTZ. Define ky^yi = k 
ii y ^ y' £ R, and ky^yi = if y ^ y' ^ R. Define q = A; if 5; — )• G i? and q = 1 

otherwise. Then the monomial corresponding to the set R' G TZs is the only monomial that 
tends to ±00 when k tends to infinity for m large enough. □ 

Example 7.7. We consider the extension of the futile cycle to incorporate two modification 
sites instead of one. The network consists of the reactions 

51 + 5*3 < ^ 55 ^ Si + 5*4 < ^ 5*7 5"! + Ss 

52 + Ss < ' Sg ^ S2 + »S'4 < ' Sq 3- 5*2 + 5*3 

fcg fcii 

with rate constants indicated next to each reaction. The enzyme catalyzes the modi- 
fication of 53 to S4 and subsequently to Sg via the formation of the intermediates 55, S7. 
Similarly, the enzyme 5*2 catalyzes the demodification of 5*8 to 5*4 and then to 5*3 via the 
intermediates Sg and Sq. A reduced basis of L-*- of this network is given by the vectors 
5*1 + 5*5 + 5*7, S2 + Sq + Sg, and53 + ^4 + 55 + + ^7 + S's + Sg. The extended rate function 
is obtained by substituting the components of with indices 1,2, and 3 by ci + C5 + C7, 
C2 + C6 + cg, and C3 + C4 + C5 + ce + C7 + cs + cg , respectively. 

The determinant of the extended rate function has the monomials 

kik^kik'jkgki2CiC2C^ and — A;2^4^6^7^io^i2CiC2C4. 

The two terms have different signs. Therefore, this network is not injective. Note that the 
degree of the monomials in the rate constant is s = 6. 

It is well known that this network can exhibit multistationarity for some choices of rate 
constants and total amounts (|16j). More generally, in [Tl] Feliu and Wiuf analyzed the 
occurrence of multistationarity in different smalls motifs accounting for enzyme sharing in 
protein modification, including the futile cycle and the two-site modification cycle. In their 
examples, all motifs that admit exactly one positive steady state for any total amounts are 
in fact injective. The motifs that can admit multiple positive steady states are obviously 
not injective. 

Example 7.8. Being injective is not a necessary condition for the existence of at most one 
positive steady state within each stoichiometric class. Consider the network with reactions 
ri: A + B -^C r2: C ^ A + B r^: A + E rr: B^O 

rz: D + E ^ F r^ - F ^ D + E r^: D ^ B + D rs: E-^0. 
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For Ri = {r2, r4, ry, rs}, we have a{Ri) = 1, while for R2 = {j"!, T3, rs, rg}, we have 
<y{R2) = —1- It follows that the network is not injective. However, by solving the steady- 
state equations together with the equations for the conservation laws, it is easily seen that 
there is exactly one positive steady state in each stoichiometric class. 



8. Networks with all steady states degenerate 

If det(Jc(/K)) is not identically zero then there exist s linearly independent reaction 
vectors — y'^ , . . . — y'^ such that y^, . . . ,y^ are also linearly independent. Therefore, 
if the dimension of the vector space 

y := {y G C\y is the reactant complex of some reaction y ^ y') ^ M" 

is strictly smaller than s then det( Jc(/k)) = 0. For example, if the reactant complexes of a 
network involve at most s — 1 species, then all steady states of the network are degenerate. 
The network with reactions {A ^ B,A ^ 0} satisfies dim 3^ = 1 < 2, consistent with our 
computation that det(Jc(/K)) = det(Jc(/K)) = 0. 

Let / G OdiM) and let 5/ = {5j G 5| z ^ 1} be the set of species with indices not in /. 
Note that Sj has cardinality s. We consider the projection of the network N to the set of 
species 5/, A/}, induced by the projection vr/ : — )• W on the coordinates not in /. For 
example, the projection of the futile cycle with / = {1,2,6} is the network with species 
Si = {53, 54, 5*5 } and reactions 

53 < ' 5*4 54 < ' ^ 53. 

If — )• occurs then the reaction is discarded and like-wise redundant reactions are re- 
moved. The matrix T{R)i varies over all sets of s reactions in A// as R varies. Similarly, 
y{R)i varies over all sets of s reactant complexes in Mi as R varies. Therefore, the re- 
quirement that det(3^(i?)/) det(r(ii)/) 7^ for some R is equivalent to the existence of 
s independent reactions in Mi such that the corresponding reactant vectors also are in- 
dependent. Since Mi has s species, a necessary condition is that Mi is open, that is, 
the stoichiometric space has dimension s. In the example above, the set of reactions 
{S3 — )• /S'5,5'5 — )• 5*4, 5*4 0} are independent and so are the reactant complexes. This 



implies (as also shown in Example 4.8) that the steady states of the futile cycle are not 
degenerate. 

The following corollary is a consequence of Corollary |7.3[ Corollary |4. 7| and the preceding 
discussion. 

Corollary 8.1. Let M = {S,C,TZ) be a network. The following statements are equivalent: 

(i) ker( Jc(/k)) n F / {0} for all k G and c G M!]:. 

(ii) det{y{R)i) det(F(i?)/) = for all sets R of s reactions in TZ and I G OdiM). 

(Hi) For all I G Od{M) and for any set of s reactions y^ — )• y'^,...,y^ — )• y''^ in the 
projected network Mi, if the vectors y^ — y'^ , . . . , y'^ — y'^ are linearly independent then 
the complexes y^, . . . ,y^ are linearly dependent. 

If any of these hold then M is not injective and any steady state is degenerate. 

If M is open, then d = and the only projection to consider is the identity. Therefore, 
condition (iii) reduces to the condition of the following corollary. 
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Corollary 8.2. Let N = {S.C.ll) he an open network. Then ker( Jc(/k)) n T / {0} for all 
K G and c G if and only if for any set of n reactions — )• y'^, . . . ,y" — )• y'" such 
that the vectors — y'"*^, . . . , — y'" are linearly independent, the complexes y^, . . . , y" are 
linearly dependent. 

If is fully open then the set of species outflow reactions provides a set of independent 
reaction vectors and independent reactant complexes. Therefore, fully open networks can- 
not fulfill that ker(Jc(/K)) n T / {0} for aU k € and c G W^. However, open networks 
that are not fully open might fulfill the condition. For example, consider the network with 
reactions 

ri: Si S3 S3 Si rs : Si + S2 —!- S3 : + 25*2 — )• S3. 

The dimension of 3^ and P agree {s = 3) and the network is open. The network has no 
outflow reactions. Reactions ri, r2 are linearly dependent and thus any set of 3 independent 
reactions must contain and r^. In that case, however, the reactant complexes are linearly 
dependent. It follows that Corollary |8.2[ ii) is fulfilled and hence all steady states of the 
network are degenerate. 

Remark 8.3. If 3^ F, in which case dim(3') < s and Corollary |8.1[ ii) is fulfilled, then 
all stoichiometric classes have either none or infinitely many positive steady states: let 
a £ be a steady state. Then there is a vector 7 G F such that 7 • y = for all 
y G 3^. Consequently, for all t G M and ^ G F"*- C 3^-*- we have e^'^^"''^^^ = 1 and hence 
= ae^'^+f is a steady state. Further, q,^ £ Gt = {c £ R!;:|log(c) - log(ae*T) G F-^} 
which intersects each stoichiometric class in exactly one point for each t [9]. It is easy to 
prove that Gt H Gf = if t 7^ t' and hence there are infinitely many positive steady states 
in each stoichiometric class. 

Remark 8.4. In ^7j, a related determinant criterion is given to decide whether or not a 
closed network can admit degenerate steady states. Fix a determinant function detr on the 



stoichiometric space F. For each choice of scalar product *c in equation (4.2) (that is, for 
each choice of c G M" ) non-singularity of Jc(/k) restricted to F is related to a polynomial 
expansion in k of detr( Jc(/k))- Each term in the determinant expansion takes the form in 



Corollary 7.6 with 'det' replaced by detr and y G y{R)i replaced by Trdy), where tTc is the 
projection onto F as defined by the scalar product *c. The criterion requires the coefficients 
of the terms in the determinant expansion to be of the same sign or zero for each choice 
of c G [3 Prop. 10.3]. It is not obvious how to check whether this criterion is fulfilled 



using computational algebra software, in contrast to the criterion in Proposition 7.6 

Remark 8.5. For fully open networks, an algorithm is provided Joshi and Shiu in |15] to 
simplify the search for coefficients of the determinant of the Jacobian that have the "wrong" 
sign. The algorithm can be applied in the present setting, that is, to closed networks, to 
restrict the sets of s reactions to consider. 

Remark 8.6. The projected networks are embedded networks as defined by Joshi and Shiu 
in jl4] . We have shown that for an injective network all embedded networks obtained by 
selecting sets of s species as above are either injective or have all steady states degenerate. 
Consequently, if a network is injective then the embedded networks with s species do not 
have the capacity to admit multiple non-degenerate steady states. 
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9. Open and closed networks and injectivity 

In [5], Craciun and Feinberg preclude multistationarity in closed networks provided that 
(1) the fully open network is injective and (2) the closed network does not have degenerate 
steady states. In a later paper [71 Th. 8.2], the authors provide a sufficient condition 
(namely that the entrapped- species projection is a normal reaction network) for (2) to hold 
provided (1) holds as well. All weakly reversible networks fulfill this condition. 

Using the results of the previous sections we now relate injectivity of a network M and 
injectivity of the associated fully open network J\f°. 

Theorem 9.1. LetM = {S,C,TZ) be a closed network. If N° is injective then the following 
statements are equivalent: 

(i) M is injective. 

(ii) As a polynomial in k, det(Jc{gK°)) has at least one monomial of total degree s. 
(Hi) ker( Jc(/k)) n T = {0} for some fixed k E and c G M+. 

(iv) There exists I G Od{M) and a set of s reactions — )• y'^, . . . , — )• y'** in the projected 
network Mj such that the two sets of vectors — y'^ , . . . ,y'' — y'^ and y^ , . . . are 
both linearly independent. 

Proof. Let {w^, . . . , o;'^} be a reduced basis of F-*- and the associated extended rate 
function for k € M^. Since N° is injective, then Jc{gK°) is non-singular for all k G and 



for any c G (Proposition 6.2). Equivalently, det( Jc(5k°)) is a polynomial in k, with all 
non-zero coefficients having the same sign and thus all non-zero coefficients of det( Jc(/k)) 
have the same sign. It follows that either (1) det(Jc(/K)) = as a polynomial in k or (2) 
det( Jc(/k)) ^ for all values of /t G and c G . The equivalence of the four statements 



follows from this observation. Theorem 5.6, Corollary |5.9| and Corollary |8.1[ □ 



Corollary 9.2. LetM = {S,C,Tl) be a closed network, J\f° = {S,C°,TZ") the associated 
fully open network, and g^ the species formation rate function of N° . Assume that N° 
is injective. Then N is injective if and only if as a polynomial in k = {ky^yi)y^yi iz-jio , 
det(Ji((7K)) has at least one monomial in kyi^yii, . . . ,kyn^yin with s reactions in IZ. If 
this is the case, then all steady states of N are non-degenerate. 

Example 9.3. Consider the open network N with reactions ri: A -\- B ^ D, r2. A ^ C, 
r3 : B ^ D. The dimension of the stoichiometric space is 3. Let ki denote the rate constant 
of reaction rj. The associated fully open network is injective, since 

det(Jc(5'K°)) = 1 + A;2 + + kiCA + kics + k2k^ + kik2CA + kik-^CB- 

However, the determinant has no monomial involving ki,k2,k3. It follows that M is not 
injective and all steady states are degenerate. Alternatively, note that y = {A, B) has 
dimension 2 < s = 3. 

In view of these results, there might exist injective networks such that their open network 
counterparts are not injective. This can only occur if some monomials with less than s rate 
constants from the true reactions have the wrong sign. We provide two examples. 

Example 9.4. Consider the network N given by the reactions: 

ri: A + C ^2A + B, r2: B + C ^ A + C, r^: A + B + C ^2B. 



PRECLUSION OF SWITCH BEHAVIOR IN REACTION NETWORKS 



23 



The stoichiometric space of M has maximal dimension 3 and thus it is open. Let ki denote 
the rate constant of reaction rj, and kA_,kB,kc the rate constants of the outflow reactions. 
Then, for k = {ki, k2, k^, kA, ks, kc), 

det{Ji{gi^)) = -2kik2ks - kik2kA - kik2kB - k2k^kB - kikAks - k^kAks + 2kik2kc 

- 2kikskc - k2kAkc + k^kAkc + kikskc - k^kskc - kAkskc- 

Since the sign of the coefficients of det(Ji(gf^)) as a polynomial in n differ, the fully open 
network, J\f°, is not injective. However, there is one monomial only in ki,k2,k3, which 
implies that M is injective and has no degenerate steady states. This network, however, 
does not have positive steady states either. 

Example 9.5. Consider the network J\f given by the reactions: 

n: A + B + C ^2A + B + 2C, r2: A + C ^ A + B + C, rs: C ^ A + B + 2C. 

The stoichiometric space of this network has dimension 2 and thus it is closed. Let ki 
be the rate constant of reaction and fix the rate constants of the outflow reactions to 
kA = ks = kc = 1. Then, for = (ki, k2, k^, 1, 1, 1), 

det(Jc(s'Ko)) = kiksCACc + kik2c\cc + kik2CACc + ks + kiCACs + kicscc - 1. 
We see that N° is not injective but since all monomials in ki,k2, k^ of total degree 2 have 



the same sign, it follows from corollaries 7.2, 7.3 and 7.6 that M is injective and there are 
no degenerate steady states. 

Remark 9.6. The open network given in [SJ §6] to illustrate that degenerate steady states 
can occur is not injective. Therefore the results of this work do not apply. 

Remark 9.7. This remark is for the readers familiar with [7J and the terminology intro- 
duced in that paper. Let be a network such that M° is injective. We have provided a 
sufficient and necessary condition for M to be injective as well, namely, that there exists 
/ G Od{M), and a set of s reactions — )• y'^, . . . , — )• y''' in the projected network A/} 
such that the set of vectors y^ — y'^ , . . . ,y^ — y'^ and y^, . . . ,y^ are both linearly indepen- 
dent. Since it is a sufficient and necessary condition, any network fulfilling the condition 
in [71 Th. 8.2] for J\f to be injective, that is, the entrapped-species projection is a normal 
reaction network, also fulfills our condition. In particular, weakly reversible networks are 
normal and hence they also fulfill our condition. 

10. Weakly Sign Determined (WSD) networks 

A square matrix is said to be a P-matrix if all principle minors are positive. If the 
principle minors are non- negative the matrix is said to be a Po-^iatrix. In [2], a criterion 
(namely that the stoichiometric matrix is WSD) is introduced that relates to injectivity 
of a network M. They restrict the class of networks to non-autocatalytic networks (NAC), 
which implies that the same species cannot both be a reactant and a product in the same 
reaction. We impose the same constraint in this section. 

In our notation the criterion states that (*) det{y(R)i) det{T{R)j) > for any set / of 
n — k distinct indices in {1, ... , n} and any set R = {y^ — )• y'^, . . . , y'^ — )• y"'} of k reactions 
from J\f if and only if —Jcifn) (minus the Jacobian) is a PQ-xaaX,x\yi [21 Th. 4.1, Th. 4.3]. 
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If J\f is fully open then (*) (with inflow and outflow reactions excluded from the sets R) is 
equivalent to —Jcifn) being a P-matrix O Cor. 4.2, Th. 4.4]. Using the results of [12], it 
follows that is injective. 

If M is closed and injective then Jcifn) is not a P-matrix as the rank is less than 
n. Our criterion for injectivity states that if there exists a G {0, 1} such that (**) 
(— l)'^ det(3^(P)/) det(r(P)/) > for any set / of n — s distinct indices in {1, ... , n} and 
any set R = {y^ — >• y'^, . . . , — t- y'^} of s reactions from Af, and at least one of them 
is non-zero, then det(Jc(/K)) is non-zero and J\f is injective. Assume that (**) holds. If 
all species in J\f appear in some reactant complex then the NAC assumption guarantees 
that the diagonal entries of Jcifn) iu rows d + 1, . . . ,n are non-zero and negative. Since 
the diagonal entries in rows 1, . . . , d are all equal to one, it follows that the product of the 
diagonal entries is a term in the expansion of det(Jc(/K)) and has sign (—1)''. Consider the 
matrix J*(/k) obtained from Jcifn) by multiplying the lower s rows by minus one. The 
matrix has full rank and is a P-matrix: if N fulfills (**) then all non-zero terms in the 
determinant expansion of Jcifn) have sign (—1)*. Each principle minor can be obtained as 
a sum of terms in the determinant expansion of J*(/k) divided by the (positive) diagonal 
entries of the rows not taken in the minor. Since this sum contains at least^the non-zero di- 
agonal product term, all principle minors are positive. Consequently, J*(/re) is a P-matrix. 
Note that if has full rank (for example if it is fully open) then J*(/k) = —Jcifn)- 

If Si is a species that is only in product complexes of M then the lower s rows of Jcifn) 
are zero in the i-th position. If Si is not involved in any conservation law, then the i-th 
column of Jdfn) is zero, det( Jc(/k)) = 0, and (**) does not hold. Generally, let Si-^, . . . , Si^ 
be the species of A/" that are only in reactant complexes such that the lower s rows of Jcifn) 
are zero in entries ii, . . . , im- If (**) holds and det( Jc(/k)) / then the columns ii, . . . ,im 
of Jc(/k) are linearly independent. Since the lower s x m submatrix is identically zero, 
it follows that the upper d x m matrix has rank m (and in particular m < d). As a 
consequence, we can reorder the species in S such that S'jj , . . . , 5i,^ are the first m species 
and we are guaranteed that there exists a reduced basis of F-*- with that order. Because 
m < d, with this order the diagonal entries of Jcifn) are non-zero and we can proceed as 
above. 

To sum up, if (**) holds, then there exists an order of the species of J\f such that the 
matrix obtained from Jcifn) by changing the sign of the last s rows is a P-matrix. Using 
the results of |12j we conclude that is an injective function and hence that M is injective. 

11. Concluding remarks 

In this paper we have provided a Jacobian criterion for the characterization of injective 
networks taken with mass-action kinetics. Injective networks have the important property 
that multiple positive steady states within any stoichiometry class cannot occur for any 
choice of rate constants. Further, the existence of multiple boundary solutions of a certain 
type is also precluded. Importantly, if an injective network has a positive steady state, 
then it cannot have any other non-negative steady state. 

Since injective networks are characterized by a non-singular Jacobian of the species 
formation rate function (when restricted to the stoichiometric space), other interesting 
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Figure 1. This work in relation to previous work of Craciun and Feinberg. 



properties of this class of networks are expected. For instance, in [18j it is shown that 
their steady-state fluxes exhibit a certain degree of robustness against variation in rate 
constants. 

The main novelty of this work is that injectivity of a network can be assessed directly 
avoiding the detour to fully open networks. In Figure [T] we show how our work relates to 
previous work on injectivity by Craciun and Feinberg in [H O [7j. 

The Jacobian criterion presented here can be effectively implemented using any symbolic 
computation software like Mathematica. Suggested steps for its implementation are the 
following (using pseudo-Mathematica commands): 

• Definitions: 

n: number of species, A: stoichiometric matrix F, 

c=c[l],...,c[n]: Concentration vector, v: rate vector (with concentra- 
tions) {ky^y/Cy)y^y'(Z'JZ, 

F=A.v: species formation rate function /«. 

• Conservation laws and associated extended rate function: 

P=RowReduce[NullSpace[Transpose[A]]]: find a reduced basis of 
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ind: vector of the indices of the first entry of each vector in P, 
Ftilde: remove entries ind of F and add the entries P.c. 

• Compute the determinant of the Jacobian of Ftilde: 

J=TabIe[ D[ Ftilde[[i]], cp] ], {i,l,n},{j,l,n}], 
D=Det[J]. 

• Determine the signs of the coefficients of D: 

m = MonomialList [D] : monomials of D, 

coefrs=DeIeteDupIicates[m/.{k— )'l,c— ^l}]: coefficient of each of the 
monomials (in the rate constants and c), and delete duplicates, 
Pos = Select [coefFs, Positive], Neg=Select[coeffs, Negative]: se- 
lect the positive and the negative coefficients, 
If Length[Pos]>0 and Length[Neg]=0: M is injective. 
If Length[Pos]=0 and Length[Neg]>0: M is injective. 
Otherwise: J\f is not injective. 

In our experience, this procedure works fast and reliably for not-so-big networks (at least 
up to around 15-20 species). For bigger networks, the computational cost in finding the 
determinant might be too high for a standard computer. In that case, one might construct 
the lists Pos and Neg above by computing one by one the signs of det(3^(i?)/) dct(r(i?)/) 
for all sets R oi s reactions and indices / G Oci{M). If two conflicting (different) signs are 
found, the algorithm should stop and the network is not injective. In Mathematica, the 
command Subsets[list,{s}] provides the subsets of a list containing exactly s elements. 
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